Add Hyperbolic Sainte-Marie Equations 1D by MarcoArtiano · Pull Request #135 · trixi-framework/TrixiShallowWater.jl

MarcoArtiano · 2026-03-29T10:15:38Z

The entropy conserving and well-balanced scheme is derived in https://arxiv.org/abs/2603.18978.

github-actions

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

JuliaFormatter

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Lines 221 to 222 in 470f555

    
           f2_fluxdiff = alpha_1 * equations.gravity * h_avg * b_jump + alpha_2 * 2 * p_avg * b_jump 
        
           f2_pointwise = (1 - alpha_1) * equations.gravity * h_ll * b_jump + (1 - alpha_2) * 2 * p_ll * b_jump

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Lines 226 to 227 in 470f555

    
           f4_fluxdiff = alpha_3 * equations.celerity^2 * h_avg * v_jump - 2 * alpha_2 * equations.celerity^2 * v_avg * b_jump 
        
           f4_pointwise = (1 - alpha_3) * equations.celerity^2 * h_ll * v_jump - 2 * (1 - alpha_2) * equations.celerity^2 * v_ll * b_jump

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Lines 231 to 234 in 470f555

    
           f2, 
        
           f3, 
        
           f4, 
        
           zero(eltype(u_ll)))

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Line 239 in 470f555

equations::HyperbolicSainteMarieEquations1D)

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Line 255 in 470f555

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Line 322 in 470f555

    
           e = (h_v^2 + h_w^2) / (2 * h) + h_p^2 / (2 * equations.celerity^2 * h) + 0.5f0 * equations.gravity * h^2 + equations.gravity * h * b

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Line 326 in 470f555

    
           @inline function Trixi.lake_at_rest_error(u, equations::HyperbolicSainteMarieEquations1D)

[JuliaFormatter] _{reported by reviewdog 🐶}

TrixiShallowWater.jl/src/equations/hyperbolic_sainte_marie_1d.jl

Line 332 in 470f555

end # @muladd

codecov · 2026-03-29T10:20:01Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.30%. Comparing base (a9ae27b) to head (1b9ceff).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #135      +/-   ##
==========================================
+ Coverage   99.28%   99.30%   +0.02%     
==========================================
  Files          96       98       +2     
  Lines        5146     5290     +144     
==========================================
+ Hits         5109     5253     +144     
  Misses         37       37

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

ranocha

Thanks a lot!

test/test_tree_1d.jl

ranocha · 2026-03-29T14:17:00Z

src/equations/hyperbolic_sainte_marie_1d.jl

+@doc raw"""
+    HyperbolicSainteMarieEquations1D(; bathymetry_type = bathymetry_mild_slope,
+                                       gravity,
+                                       eta0 = 0.0,
+                                       h0,
+                                       alpha = 3.0)
+
+Hyperbolic approximation of the Sainte-Marie system
+[`SainteMarieEquations1D`](@ref) in one spatial dimension


Please fix the docstring so that it matches the code here and not DispersiveShallowWater.jl.

src/equations/hyperbolic_sainte_marie_1d.jl

ranocha · 2026-03-29T14:19:52Z

src/equations/hyperbolic_sainte_marie_1d.jl

+function HyperbolicSainteMarieEquations1D(; gravity, H0 = zero(gravity),
+                                          b0 = one(gravity), alpha = 3,
+                                          threshold_limiter = nothing)
+    T = promote_type(typeof(gravity), typeof(H0), typeof(b0), typeof(alpha))
+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    celerity = alpha * sqrt(gravity * b0)


I think we should not use b0 to determine the hyperbolization parameter, since Escalante et al. use a reference water height h_0 (and the exact value of the bathymetry b can be changed arbitrarily by moving everything in the vertical direction).

H0 should also not be used since it is a reference for the surface elevation h + b, correct? Maybe we need to introduce another argument, h0, although this can be confusing...

Would it make sense to use H0 - b0 for the reference water height?

Yes - but b0 is not used otherwise. Thus, it would be more straightforward to use another parameter name for it, wouldn't it?

Ah, yes if b0 is not used otherwise I agree it makes more sense to use another parameter name.

What about using something like h_ref or eta0?

I am fine with h_ref and h_0. The docstring is talking about h_0. So why not that?

I am in favor of h_ref because with h_0 I would expect the relation H_0 = b + h_0 to hold.

ranocha · 2026-03-29T14:24:01Z

src/equations/hyperbolic_sainte_marie_1d.jl

+"""
+	flux_conservative_ec(u_ll, u_rr, normal_direction::AbstractVector, equations::HyperbolicSainteMarieEquations1D)
+
+Total energy conserving and well-balanced two-point flux by
+-  Marco Artiano, Hendrik Ranocha (2026)
+   On Affordable High-Order Entropy-Conservative/Stable and 
+   Well-Balanced Methods for Nonconservative Hyperbolic Systems
+   [DOI: 10.48550/arXiv.2603.18978](https://arxiv.org/abs/2603.18978)
+"""
+struct flux_conservative_ec{RealT <: Real}


Structs should have capital names like FluxConservativeEC. Please also adapt the docstring to describe what the parameters alpha1, alpha2, alpha3 mean, and how the flux can be applied.
It would also be possible to just pick one set of parameters that performs well and stick with it instead of implementing the whole family. I don't have a strong opinion on this. Since performance is not critical right now, it is definitely nice to have the flexibility - but we may recommend some reasonable default value.

I've fixed the structs' name and adapted the docstrings. We may recommend alpha_1 = 1.0, alpha_2 = 0, alpha_3 = 0, which corresponds to the split form of the shallow water terms and strong form of the remaining terms. Eventually also alpha_1 = 1.0, alpha_2 = 1.0, alpha_3 = 1.0 works (as informative example), which results in a split form for all the terms.

ranocha · 2026-03-29T14:25:09Z

src/equations/hyperbolic_sainte_marie_1d.jl

+# Convert conservative variables to entropy
+# Note, only the first two are the entropy variables, the third entry still
+# just carries the bottom topography values for convenience
+@inline function Trixi.cons2entropy(u, equations::HyperbolicSainteMarieEquations1D)
+    h, h_v, h_w, h_p, b = u
+
+    v = h_v / h
+    w = h_w / h
+    p = h_p / h
+
+    w1 = equations.gravity * (h + b) - 0.5f0 * (v^2 + w^2 + p^2 / equations.celerity^2)
+    w2 = v
+    w3 = w
+    w4 = p / equations.celerity^2
+    return SVector(w1, w2, w3, w4, zero(eltype(u)))
+end


The comment about the bottom topography and the code do not match.

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

…work/TrixiShallowWater.jl into ma/hyperbolized_sainte_marie

src/equations/hyperbolic_sainte_marie_1d.jl

JoshuaLampert

Looks mostly good to me. It might make sense to also implement a flux function computing the analytical flux because it enables using, e.g., LLF.

JoshuaLampert · 2026-03-30T17:23:09Z

src/equations/hyperbolic_sainte_marie_1d.jl

+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    celerity = alpha * sqrt(gravity * b0)


Since we only need the squared celerity it probably makes sense to also just store the square.

JoshuaLampert · 2026-03-30T19:55:42Z

src/equations/hyperbolic_sainte_marie_1d.jl

+function HyperbolicSainteMarieEquations1D(; gravity, H0 = zero(gravity),
+                                          b0 = one(gravity), alpha = 3,
+                                          threshold_limiter = nothing)
+    T = promote_type(typeof(gravity), typeof(H0), typeof(b0), typeof(alpha))
+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    celerity = alpha * sqrt(gravity * b0)


I am fine with h_ref and h_0. The docstring is talking about h_0. So why not that?

patrickersing

Thanks for adding this new system! That is a really nice addition to the package.
The implementation looks quite good already, I mainly had suggestions about naming conventions and comments.

Could you please also:

Add a convergence test and document convergence result in the PR.
Add equation specific unit tests such as https://github.com/trixi-framework/TrixiShallowWater.jl/blob/main/test/test_unit.jl#L26.
Create an entry to NEWS.md about the addition.

patrickersing · 2026-03-31T07:08:30Z

examples/tree_1d_dgsem/elixir_hyperbolic_sainte_marie_ec.jl

+# semidiscretization of the shallow water equations with a discontinuous
+# bottom topography function for a fully wet configuration


Suggested change

# semidiscretization of the shallow water equations with a discontinuous

# bottom topography function for a fully wet configuration

# semidiscretization of the hyperbolic sainte-marie equations for a smooth and periodic initial condition to test entropy conservation

patrickersing · 2026-03-31T07:11:13Z

examples/tree_1d_dgsem/elixir_hyperbolic_sainte_marie_ec.jl

+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.1)


The stepsize callback is defined but not included in the CallbackSet.

patrickersing · 2026-03-31T07:15:51Z

src/equations/hyperbolic_sainte_marie_1d.jl

+The additional quantity ``H_0`` is also available to store a reference value for the total water height that
+is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
+Escalante, Dumbser and Castro (2019) choose the hyperbolization parameter as ``c = \alpha \sqrt{g h_0}`` for some background water height ``h_0``.
+Thus, the hyperbolization parameter ``c^2`` is set by the keyword arguments `alpha` (``\alpha``), `gravity` (``g``), and `h0` (``h_0``).


Just a reminder to make the variable name h_0 consistent to the function argument once we settled on a name.

patrickersing · 2026-03-31T08:05:51Z

src/equations/hyperbolic_sainte_marie_1d.jl

+function HyperbolicSainteMarieEquations1D(; gravity, H0 = zero(gravity),
+                                          b0 = one(gravity), alpha = 3,
+                                          threshold_limiter = nothing)
+    T = promote_type(typeof(gravity), typeof(H0), typeof(b0), typeof(alpha))
+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    celerity = alpha * sqrt(gravity * b0)


I am in favor of h_ref because with h_0 I would expect the relation H_0 = b + h_0 to hold.

patrickersing · 2026-03-31T08:09:19Z

src/equations/hyperbolic_sainte_marie_1d.jl

+    gravity::RealT
+    H0::RealT
+    celerity::RealT
+    threshold_limiter::RealT


Since there is not wet/dry functionality implemented, you can remove threshold_limiter for the equation.

patrickersing · 2026-03-31T08:18:07Z

src/equations/hyperbolic_sainte_marie_1d.jl

+@inline function source_term_hyperbolic_sainte_marie(u, x, t,
+                                                     equations::HyperbolicSainteMarieEquations1D)


Please add a docstring to explain the usage of this source term.
This is not optional, but a fundamental part of the equations right?

patrickersing · 2026-03-31T08:27:13Z

src/equations/hyperbolic_sainte_marie_1d.jl

+end
+
+"""
+	FluxConservativeEC(alpha_1, alpha_2, alpha_3)(u_ll, u_rr, normal_direction::AbstractVector, equations::HyperbolicSainteMarieEquations1D)


The function name FluxConservativeEC is kept very general. I would prefer to stick with the standard naming convention and name it FluxArtianoEtal and FluxNonconservativeArtianoEtal.

We also use alpha when we determine the hyperbolization parameter. What do you think about changing one of the parameter names to beta to differ between them?

Well, alpha_i was used in the paper. Thus, we would at least have to describe this difference from the paper.
Moreover, alpha is already used in several places, e.g., for the shock-capturing parameter. However, I do not have a very strong opinion on this; it's your decision what you would like to have in TrixiShallowWater.jl.

patrickersing · 2026-03-31T08:33:31Z

src/equations/hyperbolic_sainte_marie_1d.jl

+In the formal limit ``c^2 \to \infty``, the hyperbolic approximation recovers the original Sainte-Marie system.
+The additional quantity ``H_0`` is also available to store a reference value for the total water height that
+is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
+Escalante, Dumbser and Castro (2019) choose the hyperbolization parameter as ``c = \alpha \sqrt{g h_0}`` for some background water height ``h_0``.


I would suggest to name this c_ref or c_hyp as we use the variable c already, when we compute the wave celerity in max_abs_speeds.

c_ref sounds good and would be in accordance with h_ref 👍

patrickersing · 2026-03-31T08:41:10Z

src/equations/hyperbolic_sainte_marie_1d.jl

+end
+
+# Helper function to extract the velocity vector from the conservative variables
+@inline function Trixi.velocity(u, equations::HyperbolicSainteMarieEquations1D)


Should this also extract the vertical velocity?

I would not do so. Right now, the velocity dimension is in accordance with the spatial dimensions of the equations. Moreover, the equations do not really model transport in the vertical direction (just resulting vertical velocity coming from the background incompressible Euler system).

patrickersing · 2026-03-31T08:44:11Z

src/equations/hyperbolic_sainte_marie_1d.jl

+    H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
+
+    return abs(H0_wet_dry - (h + b))


Suggested change

H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)

return abs(H0_wet_dry - (h + b))

return abs(equations.H0 - (h + b))

If we don't consider wetting and drying this can be simplified.

add hypsm

470f555

github-actions bot reviewed Mar 29, 2026

View reviewed changes

fmt

fd849ec

MarcoArtiano added 8 commits March 29, 2026 12:46

add tests

ab09b24

add docs for the fluxes

78326d9

fixup

c780f8e

fix type tests

1512512

fix type tests

d18e1e0

fix

02adfd7

fix

9943b3e

fis

8dd05fc

MarcoArtiano marked this pull request as ready for review March 29, 2026 12:40

MarcoArtiano requested a review from patrickersing March 29, 2026 12:41

MarcoArtiano added 2 commits March 29, 2026 15:44

change name of flux

6d2249c

fix exports

63c6de4

ranocha requested changes Mar 29, 2026

View reviewed changes

MarcoArtiano and others added 4 commits March 29, 2026 18:11

Apply suggestions from code review

ffe33c2

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

fix docs

a9fb37b

Merge branch 'ma/hyperbolized_sainte_marie' of github.com:trixi-frame…

3d0520b

…work/TrixiShallowWater.jl into ma/hyperbolized_sainte_marie

change name convenction for the struct

a835f87

github-actions bot reviewed Mar 29, 2026

View reviewed changes

src/equations/hyperbolic_sainte_marie_1d.jl Outdated Show resolved Hide resolved

src/equations/hyperbolic_sainte_marie_1d.jl Outdated Show resolved Hide resolved

MarcoArtiano added 2 commits March 29, 2026 18:46

fix tests

ac4c35d

fix docs

1b9ceff

JoshuaLampert requested changes Mar 30, 2026

View reviewed changes

patrickersing requested changes Mar 31, 2026

View reviewed changes

	f2_fluxdiff = alpha_1 * equations.gravity * h_avg * b_jump + alpha_2 * 2 * p_avg * b_jump
	f2_pointwise = (1 - alpha_1) * equations.gravity * h_ll * b_jump + (1 - alpha_2) * 2 * p_ll * b_jump

	f4_fluxdiff = alpha_3 * equations.celerity^2 * h_avg * v_jump - 2 * alpha_2 * equations.celerity^2 * v_avg * b_jump
	f4_pointwise = (1 - alpha_3) * equations.celerity^2 * h_ll * v_jump - 2 * (1 - alpha_2) * equations.celerity^2 * v_ll * b_jump

		# semidiscretization of the shallow water equations with a discontinuous
		# bottom topography function for a fully wet configuration

	# semidiscretization of the shallow water equations with a discontinuous
	# bottom topography function for a fully wet configuration
	# semidiscretization of the hyperbolic sainte-marie equations for a smooth and periodic initial condition to test entropy conservation


		alive_callback = AliveCallback(analysis_interval = analysis_interval)

		stepsize_callback = StepsizeCallback(cfl = 0.1)

		@inline function source_term_hyperbolic_sainte_marie(u, x, t,
		equations::HyperbolicSainteMarieEquations1D)

		H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)

		return abs(H0_wet_dry - (h + b))

Conversation

MarcoArtiano commented Mar 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Mar 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

ranocha left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MarcoArtiano Mar 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MarcoArtiano Mar 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

JoshuaLampert left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

patrickersing left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MarcoArtiano commented Mar 29, 2026 •

edited

Loading

codecov bot commented Mar 29, 2026 •

edited

Loading

MarcoArtiano Mar 29, 2026 •

edited

Loading

MarcoArtiano Mar 29, 2026 •

edited

Loading